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INTRODUCTION 


Crew  and  aircraft  losses  attributed  to  loss  of  control  at  high  angles  of  attack  con¬ 
tinue  to  occur  and  to  be  of  concern  to  both  the  Navy  and  Air  Force,  and  continue  to  be 
identified  as  an  area  in  which  further  research  is  needed.  The  need  exists  not  only  to 
reduce  such  losses  of  current  aircraft  but  to  minimize  those  of  future  aircraft  now  in 
development. 

In  1979,  the  Navy  appointed  an  independent  Executive  Review  Group  to  assess  the 
problem  relative  to  the  F-14.  Among  its  conclusions  the  group  identified  causal  factors 
associated  with  airframe  aerodynamics,  control  system  mechanization,  engine  asym¬ 
metries,  piloting  technique,  and  physiological  factors.  The  group  also  recommended 
high  fidelity  simulation  as  a  low  risk  means  of  addressing  this  problem. 

Realizing  the  potential  of  the  NADC  centriftige  (djmamic  flight  simulator)  for  this 
purpose,  the  Aircraft  and  Crew  ^sterns  Technology  Directorate  has  begun  a  program 
to  substantially  improve  its  simulation  capability  to  meet  the  requirements  of  a  high 
fidelity  simulator  capable  of  modeling  a  variety  of  aircraft.  This  capability  will  encom¬ 
pass  most  conventional  simulation  techniques,  and  it  will  have  the  capability  and  flexi¬ 
bility  to  model  aircraft  dynamics  in  all  flight  conditions  including  fully  developed  spins. 

As  part  of  that  program  the  Flight  Dynamics  Branch  (Code  6053)  was  asked  to 
provide  equations  of  motion  and  F-14  aerodynamic  and  control  system  data  adequate 
to  the  task.  This  report  is  intended  to  document  those  recommendations  and  to  serve 
as  a  basis  for  subsequent  simulation  software  revisions.  Since  the  intent  is  to  use  the 
simulation  as  a  research  tool,  considerable  emphasis  is  given  to  flexibility. 

It  is  considered  beyond  the  scope  of  this  document  to  implement,  in  the  form  of 
programs  and  subroutines,  the  recommendations  presented.  In  the  interest  of  com¬ 
pleteness  however,  certain  programming  approaches  are  suggested.  While  efforts  are 
made  to  make  those  suggestions  readily  programmable,  it  is  realized  that  the  sug¬ 
gested  methods  may  not  necessarily  be  the  most  computationally  efficient.  While  some 
changes  may  be  required,  it  is  felt  that  implementation  of  all  major  features  can  be 
preserved  in  the  final  software  implementation. 
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EQUATIONS  OF  MOTION 


NONLINEARIZED  DYNAMIC  EQUATIONS 

It  is  the  need  for  the  proposed  simulation  to  represent  the  extreme  high  angle  of 
attack  flight  regime  which  distinguishes  it  from  other  simulations.  Therefore  the 
requirements  of  this  feature  are  the  driving  factors  in  settling  on  the  equations  of 
motion.  Several  considerations  peculiar  to  this  type  of  simulation  may  be  listed  as: 

1.  Nonlinearity  of  data 

2.  Uncertainty  about  form  of  data 

3.  Importance  of  various  asymmetries 

4.  Large  angle  motions 

Just  how  each  of  these  factors  impacts  the  form  of  the  equations  of  motion  will  be 
pointed  out  in  this  section. 

The  aerodynamic  force  and  moment  data  typically  changes  very  rapidly  and  often 
unexpectedly  with  changes  in  the  orientation  of  the  free  stream  velocity  vector.  Fur¬ 
thermore,  there  are  no  guarantees  that  trim  is  possible  at  all  flight  conditions  of 
interest. 


The  aircraft  may  undergo  such  extremely  violent  excursions  that  any  state  of 
equilibrium  may  be  simply  a  very  brief  transient  state.  Under  these  conditions  the 
conventional  stability  derivatives  have  little  meaning.  That  is,  the  concept  of  air¬ 
craft  motions  consisting  of  small  perturbations  from  an  equilibrium  condition  is  not 
applicable.  For  this  reason  the  usual  linearization  of  the  rigid  body  dynamic  equations, 
whether  for  analysis  or  simulation,  is  inappropriate.  The  full  nonlinear  rigid  body 
dynamic  equations  must  be  used. 


These  equations  arise  from  the  six  degree  of  freedom  equations  of  Newtonian 
dynamics. 


T  =  I  n 
dt  n 


which  define  force  and  moment  as  the  first  derivative  of  linear  and  angular  momentum, 
respectively.  These  equations  are  defined  for  inertial  (nonaccelerating,  nonrotating) 
reference  frames. 
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If  we  wish  to  write  the  force  and  moment  equations  in  a  reference  frame  that  may 
be  rotating,  we  must  remember  that  di^erentiation  of  any  vector  V  in  such  a  system 
requires  an  additional  term. 


dV  ^  _i 

-:jr=v+w  ^xv 

dt  co-ord 


to  represent  the  inertial  quantity. 


The  assumption  that  mass  and  inertia  are  approximately  constant  for  the  time 
required  to  pass  through  the  equations  of  motion  once  yields. 


F  =  ma* 


T  =  I  n 
n 


+  (nx  l^fi) 


Making  this  assumption  essentially  selects  body  axes  as  the  reference  frame  for  all 
forces  and  moments  because  it  is  the  only  reference  frame  in  which  inertias  are  con¬ 
stant  with  time.  To  avoid  unnecessary  axes  conversions,  it  will  be  assumed  that  all 
aerodynamic  data  will  be  in  body  axes. 


If  the  aircraft  motion  is  considered  to  be  a  combination  of  translation  of  the  center 
of  mass  and  rotation  about  the  center  of  mass  then  all  forces  may  be  considered  to  act 
through,  all  mass  may  be  considered  concentrated  in,  and  the  origin  of  the  body  axis 
system  may  be  attached  to,  the  center  of  gravity. 

By,  differentiating  the  expression  for  position  in  body  axis  twice,  using  the  pre¬ 
scribed  differentiation  formula,  we  may  arrive  at  the  standard  dynamic  equation  for 
inertial  acceleration. 


V  n  X  V 


Since  the  center  of  mass  is  always  at  the  origin,  R  =  0,  and  body  axis  velocity  is  true 
inertial  velocity,  its  components  are  defined  as 

-A 

=  ui  +vj  +wk 

But  body  axis  acceleration  is  related  to  inertial  acceleration  by 

-X  ^  ^  ^ 

a  =  V  +  (  fix  V). 

*  a  is  defined  as  inertial  acceleration. 
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Similar  arguments  apply  to  show  that  the  body  axis  components  of  angiilar  velocity  are 
the  components  of  the  true  inertial  rate. 

i  a  jk 

n  =  pi  +  qj  +  rk. 

Angular  acceleration  ^in  body  axis  is  also  the  inertial  angular  acceleration  since  the 
coriolus  terms  ( fix  il)  are  zero  by  the  properties  of  the  vector  cross  product. 

.X  ^  A  A 

“  =  pi  +  4i  +  rk. 

Expanding  the  linear  acceleration  equation  gives  the  more  common  form 

u  =  a^  +  (rv  -  qw) 

V  =  +  (pw  -  ru) 

w  =  a^  +  (qu  -  pv) 

Because  accelerations  can  be  easily  integrated  to  obtain  rates  and  positions,  it  is 
most  useful  to  write  the  force  and  moment  equations  in  terms  of  accelerations. 


It  is  still  a  matter  of  choice  as  to  which  forces  and  moments  are  included  in  the 
equations.  The  significant  forces  may  be  categorized  as  aerodynamic,  thrust,  and 
gravitational,  making  the  rectilinear  equations 


=  g/W^QSZ  C  . 

+  C0S€ZT. 

1.  XI 

1. 

=  g/w  |Qsr 

1  +  g  cos  9 

=  g/w{QSZ  C  . 

t  Z1 

+  sin€  ZT. 

1 

sin  9 


cos  9  COS0. 


Similarly  the  moment  equation 


T  =  I  n  +  (n  X  I  n ) 
n  '  n  ' 


may  be  written 
LJ 

n 


or 


^  “  '‘12  ^'x  *“2  "y  *'*23  “z 
'  ■  '*13  “x  *<*23  ‘"v  *'*3  “z 
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Where  Mx,  My,  include  aerodynamic  moments,  thrust  moments,  corrections  for 
center  of  gravity  shift,  as  well  as  inertia  coupling  moments,  and  engine  gyroscopic 
moments. 

If  the  applied  force  components  are  defined  as: 

F  =  Qs  2  C  .  +  cos  e  I  T. 

X  XI  1 

F  =  Qs  2  C  . 

y  yi 

F  =  Qs  2  C  .  +sine  2T. 
z  zi  1 

the  form  of  the  moment  terms  are  kept  more  manageable. 

M  =  QSb  2  C,.  +  (I  -  I  )  qr  +p  (ql  -  rl  ) 

X  li  y  z'  ^  xz  xy' 

2  2 

+  (q  -  r  )  I  +  sin  c  z  1  .  T. 

yz  yi  1 


+  F  (Az  J  -  F  (Ay  ) 
y  '  cg^  z  '  ^cg' 

+  q  sin  e  Z  I  w  . 

e  ei 


M  =  QS  c  2  C  .  +  (I  -  I  )  pr  +  q  (rl  -  pi  ) 

y  mi  2  x'  ^  ^  '  xy  yz' 

+  (r^  -  p^)  I  +2  1.  T. 

'  '  xz  Ti  1 


+  F  (Ax  )-  Fx  (Az  J 
z  '  eg'  eg' 


-  r  cos  e  2  I  w  .  -  p  sin  e  Z  I  «  . 

e  ei  e  ei 


M  =  QSb  Z  C  .  +  (I  -  I)  pq  +  r  (pi  -  ql  ) 
z  m  X  jr  yz  xz' 

2  2 

+  (p  -  q  )  I  +  cos  e2  1  .  T. 

xy  yi  1 


+  F  (Ay  )-  F  (Ax  J  +  q  cos  e  Z  I  ^  . 
X  '  “^cg  y  eg'  e  ei 
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MASS  AND  INERTIA  ASYMMETRIES 

An  important  factor  in  the  choice  of  these  equations  is  the  ability  to  account  for 
mass  and  engine  asymmetries. 

To  this  end,  the  eg.  may  be  shifted  in  any  of  three  directions  from  its  nominal 
location  ^ycg>  ^cg^*  products  of  inertia  are  assumed  zero,  and  each 

thrust  and  engine  moment  arm  is  accounted  for  separately. 

In  order  to  account  for  the  weight  and  inertia  contribution  of  individual  aircraft  or 
store  components,  these  are  expressed  as  the  summation  of  terms  some  of  which  may 
be  functions  of  variables  such  as  switch  positions,  wing  sweep,  or  time. 


I  = 

S 

I  . 

I 

=  21  . 

X 

XI 

xy 

XZl 

I  = 

2 

I  . 

I 

=  21  . 

y 

yi 

yz 

yzi 

I  = 

2 

I  . 

I 

=  21  . 

z 

Z1 

xz 

xzi 

W  =  2  W. 

1 


The  numbers,  dj,  represent  the  entries  of  the  inverse  of  the  inertia  matrix 


Since  the  product  of  the  inertia  matrix  and  its  inverse  are  the  identity  matrix. 


We  can  write 
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and  finri  each  entry  by  Cramer's  rule  which  for  gives 


1  “  det  I 


1 

0 

0 


-I 


xy 


-I 


yz 


-I 


xz 


-I 


yz 


In  doing  so  we  obtain  the  relations 

det  I  =  I  I  I  -  2  I  I  I  -  (I  1  ^+11  ^+11 

n  X  y  z  xy  xz  yz  x  yz  y  xz  z  xy  ' 

d,  =  (I  I  -  I  ^)/det  I 

1  '  y  z  yz  n 

d„  =  (I  I  -  I  ^)/det  I 

2  '  X  z  xz  '  n 

d„  =  (I  I  -  I  ^)/det  I 

3  '  X  y  xy  '  n 
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=  (I  I  + 1  I  )/det  I 
'  yz  xz  xy  z  n 

=  (I  I  + 1  I  )/det  I 
'  xy  zy  xz  y  n 

=  (I  I  II  )/det  I 
'  xy  xz  —  yz  X  n 


LARGE  ANGLE  MOTION  CONSIDERATIONS 


Since  the  simulation  may  be  required  to  travel  through  very  large  angular  motions 
it  is  essential  to  describe  both  wind  vector  and  earth  axis  orientation  so  that  all  pos¬ 
sible  orientations  are  algebraically  definable  and  continuous. 

Velocity  vector  orientation  is  definable  any  time  total  velocity  is  not  zero.  Its 
magnitude  is  given  by: 


=  4 


2  ^  2  ^  2 
U  +  V  +  W 


If  ^  0,  sideslip  is  defined  as 


0  =  sin 


7 


NADC-80220-60 


which  is  continuous  over  the  range  -90  ^  <  90*.  Angle  of  attack  is  definable  any 

time  the  velocity  vector  has  a  projection  in  the  XZ  body  axis  plane.  That  is  if 

2  2  ,  ^ 
u  +w  p  0. 


The  normal  definition  for  angle  of  attack  however. 


u 

is  not  defined  ato  =  ±90*  and  further  can  not  distinguish  between  the  range  0*<a<90* 
and  the  range  -90*<a<180*,  nor  the  range  -90*<a<0  and  the  range  90*<a<180*. 

An  alternate  expression  valid  over  the  entire  range  -180“  <  a  <180*  (not  inclusive)  is 


a  =  2  tan 


-1 


(  w 

u 

2  2 

U  +  w 

For  quadrants  I  and  n  this  expression  follows  directly  from  the  geometry  below. 


QUADRANTS  III  &  IV 


For  quadrants  HI  and  IV  the  geometry  says 

f  /9  ^  ^U^  +  W^ 

tan  a/2  =  _ ! _ 

w 

which  can  be  altered  to  the  desired  form  by  multiplying  by 

J~2  2 

u  +  ^  u  +  w 
u  ^  u  +  w 
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V2  2 

u  +w  /0,a/±180®.  Conversely  if 


u  =  -  Vu+v:  ,  «=  180®. 

Differentiating  the  expressions  defining  angle  of  attack  and  sideslip  and  substituting 
the  relations  on  page  4  between  apparent  body  axis  and  inertial  accelerations  gives 
expressions  for  rates, 

au-aw-v(pu+  rw) 

i=— 5 - ^ -  +q 

2  2  ^ 

(u  +  w  ) 

2  2  2 
.  a  (u  +  w  )  -  V  (a  u  +  a  w)  +  V._  (pw  -  ru) 

(j=-i - ^  i _ 

Vu  + 

T  ’ 

in  terms  of  readily  available  quantities. 

Similar  problems  occur  in  representing  earth  axis  orientation  with  standard  Euler 
angles  which  are  obtained  by  integrating  the  Euler  angle  rates 

0  =  p  +  (q  sin  0  +  r  cos  0  )  tan  S 

0  =  (q  sin  0  +  r  cos  0  )/cos  8 

0  =  q  cos  0  -  r  sin0. 

At  8  =±90®  the  expressions  for  0  and  0  become  undefined  and  their  values  ambiguous. 

One  possible  way  of  avoiding  this  problem  is  by  use  of  quaternions  which  are  both 
continuous  and  unambiguous.  Since  this  method  was  previously  used  for  the  ACM 
simulation  no  elaboration  is  necessary  here. 

Should  use  of  an  alternate  method  become  necessary,  one  is  suggested  by  the  fact 
that  if  the  Euler  rates  are  known,  no  such  problem  exists. 

p  =  |J>  -  0  sin  8 

q  =  8  cos  8  +  4/  cos  8  sin  0 

r  =  0  cos  8  cos  0  -  8  s'm  <t> 

This  directional  quality  suggests  that  a  set  of  inverse  Euler  angles  be  used  to  deter¬ 
mine  the  orientation  of  earth  axes  relative  to  body  axes.  The  inverse  Euler  angles 
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are  defined  in  the  same  manner  as  their  conventional  counterparts,  except  that  the 
sequence  of  rotations  begins  at  body  axes  and  ends  with  earth  axes  rather  than  vice 
versa. 


-1 


=  r 


-1 


=  q  cos  -  p  sin  <li  ^ 


^  =  (p  cos  'i'  +  q  sin  >1/  j)  cos  0  ^  - -  ^ 


-  r  sin 


This  equation  is  not  singular  for  any  orientation.  Furthermore,  the  matrix  which  maps 
to  earth  axes  using  inverse  Euler  angles  is  identical  in  form  to  that  which  maps  earth 
axes  to  body  axes  using  the  standard  set. 


cos  ^  cos  e 
^cos  <1/  sin  8  sin  <t) ' 
-cos  0  sin  0 
'  cos  0  sin  6  cos  <t> 
+sin  0  sin  0 


sin  0  cos  e 
^sin0  Sind  sin  0' 
+COS0  cos  0 
'  sin  0  sin  8  cos  0  ’ 
-cos  0  sin  0 


-sin  8 


cos  8  sin  0 


cos  8  cos  0 


Since  each  transformation  is  between  orthogonal  coordinate  systems,  the  inverse 
of  the  above  matrix  is  its  transpose.  Because  the  above  transformation  and  the  iden¬ 
tical  transformation  with  the  inverse  Euler  angles  substituted  are  inverse  operations, 
the  two  sets  of  angles  may  be  related  by  equating  forms  across  the  diagonal  of  the 
above  matrix. 


0  =  -  sin 


-1 


0  =  cos 


-1 


cos  0  , 

sin 

8  , 

’  -1 

-1 

cos  0 

cos 

8  , 

-1 

-1 

cos  8 


0  =  cos 


cos  8  ^  cos  0  ^ 
cos  8 


To  express  the  inverse  angles  in  terms  of  the  standard  angles  simply  reverse 
their  positions  in  this  form. 
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While  the  expressions  for  the  standard  Euler  angles  is  still  undefined  at  fl  =  90®, 
the  inverse  angles  do  not  experience  this  problem. 

Bates  and  positions  may  be  calculated  by  numerical  integration.  An  easy  method 
of  implementation  is  to  model  each  integrator  as  an  element  of  a  linear  system  using 
the  digital  representations  to  be  described  later.  Those  variables  needing  integration 
are; 


II 

p  =  J  fidt 

1 

II 

> 

q  =  y'qdt 

1 

II 

r  =  J"  fdt 

If  alternate  Euler  angles  are  used  to  represent  earth  axis  orientation  we  have  also 

*-i  °  I *-i 

»-i  -  /"-i 

“  J" 

Standard  Euler  angles  may  be  determined  by  direct  conversion.  Translational 
rates  may  be  converted  to  earth  axes  using  the  matrix  on  page  12,  and  then  integrated 
to  keep  track  of  position  and  altitude 


X  = 

dt 

J 

'  e 

Y  = 

/'v 

dt 

j 

f  e 

h  = 

dt 

e 


2 

Using  a  table  look-up  to  determine  Q/m  and  speed  of  sound,  a,  from  atmospheric 
data  allows  the  calculation  of 


M  =  V^a 

Q  =  (Q/m^) 
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DATA  FORM  AND  USE 

An  important  fact  concerning  extreme  high  angle  of  attack  and  spin  simulation  is 
that  the  aerodynamic  phenomena  which  generate  the  forces  and  moments  are  not  easily 
predictable  or,  in  some  oases,  well  understood.  The  significant  phenomena  may  change 
radically  from  one  aircraft  to  another,  or  from  one  angle  of  attack  or  rotation  rate 
range  to  another.  In  addition,  what  phenomena  are  known  to  exist  on  any  given  air¬ 
craft  is  something  which  changes  with  time  as  new  test  techniques  are  developed  or  as 
new  information  is  extracted  from  flight  test  data. 

By  providing  data  storage  and  equations  of  motion,  we  are  faced  with  a  great  deal 
of  uncertainty  about  exactly  which  aerod3mamic  coefficients  are  significant,  to  what 
extent  they  may  be  linear  or  nonlinear,  and  even  of  what  variables  they  may  be  func¬ 
tions.  Essentially  what  is  needed  is  a  programming  scheme  flexible  enough  to  accept 
a  wide  variety  of  data  forms  without  reprogramming. 

A  nonlinearized  form  of  the  equations  of  motion  have  been  chosen  so  that  in  the 
absence  of  linearizing  assumptions,  either  linear  or  nonlinear  data  may  be  used.  The 
aerodynamic  data  has  further  been  resolved  into  body  axis  components  and  expressed 
in  standard  force  and  moment  coefficient  form.  Each  component  is  expressed  simply 
as  a  term  in  a  summation.  This  implies  the  generalization  of  data  look  up  process 
which  frees  the  aerodynamic  and  thrust  data  of  any  specified  form. 

Under  such  a  programming  approach  each  aerodynamic,  thrust,  and  some  weight 
and  inertia  terms,  would  be  represented  by  a  data  table  of  up  to  three  dimensions. 

The  table  would  contain  all  ordinates  and  their  corresponding  data  points,  but  it 
would  also  contain  a  number  of  control  integers.  These  integers  would  specify  a 
number  of  operations  to  be  performed  on  the  result  of  the  raw  table  look  up  function. 

Such  integers  would  perform  the  following  functions; 

1.  Specify  the  number  of  ordinates  in  the  table 

2.  Identify  each  ordir  o  of  the  table 

3.  Identify  the  equation  in  which  the  table  represents  a  term 

4.  Specify  any  variables  by  which  the  result  of  the  table  look  up  must  be  multi¬ 
plied  before  adding  to  the  equation  of  motion  (for  linear  data) . 

To  accomplish  this  it  will  probably  be  necessary  to  establish  a  coding  system  whereby 
any  variable  (such  as  accelerations,  rates,  positions,  load  factors,  mach  number 
altitude,  dynamic  pressure,  etc.)  or  equation  is  assigned  an  integer,  which  identifies 
it  to  the  program  and  allows  it  to  be  specified  by  entering  its  number  in  the  input  data 
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set.  This  coding  will  also  be  necessary  for  specification  of  input  and  output  variables 
for  the  control  system.  It  may  even  be  possible,  and  probably  preferable,  to  use 
Alpha  numeric  character  sets  rather  than  numerical  codes. 

B  is  possible  that  any  data  table  may  be  applicable  only  over  a  specified  range  of  a 
particular  variable.  For  instance,  a  given  aerodynamic  coefficient  may  be  highly 
linear  over  a  given  range  of  ordinates  and  nonlinear  outside  that  range.  Therefore,  it 
is  necessary  to  be  able  to  specify  a  number  of  variables  and  limiting  values  outside  of 
which  the  table  is  bypassed  and  "cut  off".  Such  cut-off  variables  may  not  always  be 
restricted  to  those  associated  with  the  table  ordinates. 

When  these  features  are  imolemented,  it  will  be  possible  to  add  or  delete  virtually 
any  linear  or  nonlinear  data.  The  capability  exists  to  completely  reprogram  the  form 
of  the  aerodynamic  data  simply  by  changing  data  tables.  TTiis  capability  is  essential 
to  the  flexibility  required  to  maintain  the  usefulness  of  the  simulation  software  in  an 
uncertain  and  rapidly  changing  technology  area. 

EQUATION  SUMMARY 

Linear  Degrees  of  Freedom; 

a  (Fyw-  sinfl) 

X  X 

a  =  g  (F /W  +  sinp  cos  fl) 

y  y 

a  =  g  (F  /W  +  cos  0  cose) 
z  z 

Rotational  Degrees  of  Freedom;  ‘ 


F  =QSSC  . 

y  yi 

F  =QS2C  .  +sine  ST. 
z  Zl  1 
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Body  Axis  Moment: 

M  =  QSb  2  C„ .  +  sin  Z  «  .  T. 

X  ei  ‘'yi  1 

^y  ^e^ei 

2  2 

+  (I  -1  )qr  +p(ql  -rl  )  +(q  -  r  )I 
'  y  z'  xz  xy  yz 


M  =QSc  2  C  .  +2:j_.T. 
y  mi  Ti  1 


+  F  (Ax  J  -  F  (Az  )  -  (r  cos  e+psine)2l  w  . 
z  '  cr  X  '  eg  e  ei 


2  2. 


+  (I^-y  Pr  +  q  (r  I  -  pi  )  +  (r  -  P  ) 


xz 


M  =  QSb  2  C  .  +  cos  €  2e  .  T. 
z  m  yi  1 


+  F  (Ay  J  -  F  (Ax  J  +  q  cos  1 2  I  w 

X  '  cg^  y  eg'  e  ei 

2  2 

+  (I  -I  )  pq  +  r  (pi  -ql  )  +  (p  -  q  )  I 

V  X  y'  yz  ^  xz'  '  '  Xy 


Weight  and  Inertia: 

Ix  “  “Ixi  ly  ~  “  2Izi 


Ixy'^Ixyi  lyz  ^^zi  ^xz  “^xzi 


W  =  2  W. 

1 

det  I  =  I  I  I  -  21  I  I  -  (I  I  ^+11  ^+11 

n  X  y  z  xy  xz  yz  x  yz  y  xz  z  xy  ' 


d 

d 

d 


1 

2 

3 


(I  I 

-  I 

y  z 

yz 

(I  I 

-  I  ' 

X  z 

XZ 

(I  I 

- 1 

'xy 

xy 

)/det  I 

n 

)/det  I 

n 

)/det  I 

n 
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d,  „  =  (I  I  + 1  I  )/det  I 

12  '  yz  xz  xy  z  n 

d,„  =  (I  I  +I  I  )/det  I 

13  xy  zy  xz  y  n 

d„„  =  (I  I  +I  I)/detI 

23  xy  xz  yz  x^  n 

Velocity  Vector  Orientation: 


,222 

V^  =  V  u  +  V  +  w 


IFV^  0, 


IF  V^  /  0 


-1  /  V 


0  =  sin 


2  2 

IFu+w  =  0,  a  =  a  a  = 

LAST’ 


a  ,  0=  0 

LAST’  LAST 


IF 


IF 


J2  ...2 

u=  -^u 


+w  ,  or  =  180 

/  J~2  ^ 

u  ^  -yu  +w 
-1 


a  =  2  tan 


w 


u  +  yu' 


2  2 
+  \v  ) 


a  u  -  a  w  -  V  (pu+nv) 

*  Z  X 

a  =  — '  ■■  +  Q 

2  2  ^ 
(u  +  w  ) 


NADC-80220-60 

Earth  Axis  Orientations  by  Quaternions 
Body  axis  rate  derivatives 
u  =  a^  +  (rv  -  qw) 

V  =  a^  +  (pw  -  ru) 
w  =  a  +  (qu  -  pv) 

Z 

Obtain  by  integration 
u,  V,  w,  p,  q, 

X.  Y.  h 

Atmosphere  dependent  quantities 
M  =  V^a 

Q  »  (Q/M^) 
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CONTROL  SYSTEM  DIGITIZATION 


APPROACH 

Previous  simulations  on  the  NAVAIRDEVCEN  centrifuge  employed  an  analog 
representation  of  the  aircraft  flight  control  system.  For  reasons  of  run  efficiency, 
reliability,  and  repeatability,  and  for  the  capability  to  rapidly  change  control  systems 
as  well  as  to  transfer  the  simulation  to  other  Navy  simulation  facilities,  a  digital 
representation  of  the  aircraft  flight  control  systems  is  required.  The  following 
development  outlines  a  flexible  approach  for  representing  such  systems  on 
NAVAIRDEVCEN  simulation  facilities. 

Inputs  to  the  control  system  may  be  considered  to  be  any  of  a  number  of  control 
deflections  from  various  cockpit  devices  and  any  of  a  number  of  parameters  which 
define  the  flight  condition  of  the  aircraft.  Outputs  may  be  surface  deflections,  thrust 
level,  or  stick  force  feedback  terms.  The  output  terms  may  be  used  either  for 
determining  forces  and  moments  acting  on  the  aircraft  or  for  the  operation  of  cockpit 
hardware. 

All  elements  of  the  system  may  be  classified  as  either  linear  or  nonlinear.  A 
linear  element  may  be  defined  for  this  purpose  as  one  described  by  a  transfer  function. 

Nonlinear  elements  are  those  described  by  nonlinear  functions  in  which  time  is  not 
an  explicit  variable,  i.  e. ,  nonlinear  differential  equations  are  not  considered.  In 
general,  any  control  system  may  be  viewed  as  one  or  more  networks  of  linear  ele¬ 
ments,  separated  by  one  or  more  nonlinear  elements  and  joined  together  in  a  par¬ 
ticular  element  structure. 

A  general  control  system  may  have  any  number  of  arbitrarily  connected  elements. 
Therefore,  it  would  be  extremely  difficult  to  devise  a  completely  general  method  to 
model  all  possible  control  systems. 

tt  is  possible  to  represent  networks  of  linear  elements  in  a  very  compact  and 
general  fashion  using  state  space  notation,  and  to  classify  most  nonlinear  elements 
as  one  of  a  few  types.  The  nonlinear  element  structure  is  almost  always  peculiar 
to  a  particular  control  system  and,  if  altered,  will  unavoidably  require  some  degree 
of  reprogramming.  By  altering  numerical  values  which  are  obtained  from  input  data 
tables,  significant  changes  may  be  made  without  reprogramming. 

In  addition,  many  different  nonlinear  element  structures  may  be  stored  as  sub¬ 
routines,  which  may  be  conditionally  called,  allowing  rapid  change  of  control  systems. 


If  (ICS.  EQ.  1)  CALL  CONSYS  1 


etc. 
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The  control  system  may  then  be  conveniently  partitioned  into  successive  nonlinear 
and  linear  stages. 

Ibe  boundaries  of  a  nonlinear  stage  are  established  by  encountering  either  a 
linear  element,  or  the  output  of  a  linear  element  not  contained  in  any  previous  linear 
stage,  while  proceeding  in  the  direction  of  signal  flow. 

It  is  necessary  to  preserve  several  variables  used  internally  to  the  control  system 
to  link  linear  and  nonlinear  stages.  The  control  system  will  not  generally  be  e^ressible 
as  a  single  linear  system.  Hiis  type  of  reduction  is  usually  performed  assuming  all 
nonlinear  test  functions  are  positive. 


Most  nonlinear  elements  may  be  categorized  as  one  of  three  types. 

Schedule  -  Any  nonlinear  function  requiring  data  interpolation 

or  data  look  up. 

Computational  Nonlinearity  -  A  nonlinear  function  expressible  algebraically. 

Switching  Nonlinearity  -  Any  logical  test  including  on-off  and  limiting. 

A  schedule  mjy  represent  an  element  such  as  a  nonlinear  gearing,  or  it  m^ 
represent  programmed  configuration  changes  such  as  wing  sweep,  or  it  may  be  used 
to  actively  change  control  system  parameters;  limits,  transfer  function  gains,  or 
coefficients  of  computational  nonlinearities. 

Any  linear  stage  may  be  considered  a  network  of  linear  elements  represented  by 
a  transfer  matrix  equation. 
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(s) 

Gii(s)  ••• 

(S) 

Y2(s) 

• 

II 

G21<^)  S2<">  -  Sk<") 

•  •  '  • 

U^fs) 

• 

• 

•  •  • 

•  •  • 

°J1  <=>  °J2  »’>  °Jk 

■ 

• 

Uk(s) 

Use  of  the  state  space  notation  to  calculate  the  system  outputs  minimizes  compu¬ 
tational  requirements. 

Y(t)  =  CX(t)  +  Du(t) 

X(t  +  T^)  =  A*X(t)  +  Bu(t) 

This  state  space  form  implies  both  a  digital  representation  of  the  entries  in  the 
transfer  function  matrix,  as  well  as  conversion  to  a  more  compact  form.  Appendix  A 
provides  details  of  how  to  digitally  represent  continuous  elements,  while  appendix  B 
explains  how  to  apply  state  space  notation  to  achieve  the  above  form.  The  salient 
features  of  appendices  A  and  B  will  be  summarized  here. 

The  outputs  of  each  stage  msy  be  efficiently  calculated  in  real  time  knowing  only 
the  entries  of  the  two  state  space  matrices  C  and  D  and  the  input  and  state  vectors. 

The  state  vector  m^  be  updated  between  successive  calculations  of  the  output 
vector  if  the  A*  and  B  matrices  are  known.  The  entries  of  the  state  space  matrices 
may  be  determined  directly  from  the  entries  of  the  transfer  function  matrix  by  specify¬ 
ing  the  following. 

For  each  linear  stage; 

1.  No.  of  nonzero  transfer  matrix  entries. 

2.  Order  of  highest  order  entry. 

3.  Number  of  outputs. 

4.  Number  of  inputs. 

For  each  transfer  function; 

1.  Number  of  first  order  numerator  factors. 

2.  Number  of  second  order  numerator  factors. 
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3.  Number  of  first  order  denominator  factors. 

4.  Number  of  second  order  denominator  factors. 

5.  Transfer  function  gain  (const.). 

6.  First  order  time  constants,  second  order  frequencies  and  damping  ratios. 

7.  Location  in  the  transfer  matrix  (subscripts). 

Each  transfer  function  is  first  expanded  into  two  polynomials  in  S  by  an  iterative 
procedure  which  calculates  the  coefficients  of  any  polynomial  expressed  in  factored 
form.  If  8  is  temporarily  redefined  as  the  number  of  factors  which  have  been  multi¬ 
plied  out,  the  polynomial  after  8  terms  is 


p  =  E 


1=0 


where 


a  =1 
oo 


All  complex  factors  must  appear  in  conjugate  pairs  to  assure  all  a^j  are  reaU 
Therefore,  if  the  next  factor  is  first  order 

^(8+1)1  ^  ^^8  (l-l)  ^8i 

while  noting 


a„.  =  0  if  i  <  0  or  i >  8 

8i 


If  the  next  term  is  second  order. 


a  =('a<^+2{u;a  +a  ) 

(8+2)  i  '  8i  n  '  n  e(i-l)  8  (i-2)^ 


The  expansion  process  is  continued  until  all  n  numerator  coefficients,  b  .  and  all 

denominator  factors,  a  ,,  are  obtained, 

ml 

A  matrix  of  coefficients  reflecting  the  expansion  of  the  first  order  factors 


m 

h  z'’  =  (z-l)^(z+l) 
j=o  ^ 


m-i 
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For  m  >  1  >  o  comes  about  as  a  result  of  using  the  Tustin  transformation. 

S  =  ^ 

T3  (z+l) 

To  produce  a  difference  equation 


m  P  m-1  P 

Y(t)  =  E  k  -JL  U  (t-(m-j)  T  )  -  £  _iL  Y  (t-(m-j)TJ 


J==o 


3=0  P 


2m 


2m 


Whose  coefficients 


P.. 

1] 


n 


z 


_m-i 


h 


a  2^ 
mi  s 


h.. 

13 


have  the  relationships 


a. 

1 


P 


P 


2  m-i 
2m 


b. 

1 


k 


1  m-i 
P 

2m 


With  the  state  space  matrices 


JTT  7a£. 
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Where  a*  =  a^.  is  selected  from  the  ith  row  of  the  matrix 
J  ji 


a  a  a 
11  21  ...  ml 


a  a  a 
12  22  ...  m2 


a  a  a 
18  28  ...  m8 


one  row  for  each  of  the  8  elements. 


When  multiplying  the  jth  row  of  A*  by  the  ith  column  of  the  state  matrix. 

X, 


X  = 


^11  '^12 


18 


21  22  ...  28 


X  ,  X  „  X  , 
ml  m2  . . .  m8 


Each  column  of  X  is  the  state  vector  for  each  of  the  8  elements  in  the  system. 


Also 


B'  = 


®11  ®12...®18 


®21  ®22  ...®2e 


6  H  B 

ml  m22  ...  m8 


W1\ere  for  each  of  the  8  elements  each  column  is  defined  by. 

B  =  b.  -  a.  b  m  >  i  >  1 
i  1  1  o  ~ 

The  matrices  C  and  E  are  found  such  that  If  the  8th  element  occupies  the  jth  row  and 
kth  column  In  the  transfer  matrix 

And  if  =  1  otherwise  E.^  =  o 


Co.  =  1 


c„  =  o. 
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And  the  matrices  B  and  D  are  obtained  by  matrix  multiplication 

B  =  B'  E  and  D  =  C  b  E 

o 

Where 


All  linear  stage  schematics  may  be  reduced,  through  block  diagram  algebra,  to 
a  transfer  matrix  form. 

Only  the  multiplication  of  the  state  space  matrices  must  be  accomplished  in  real 
time.  The  matrix  entries  may  be  calculated  during  program  initialization.  The  state 
matrix  should  normally  also  be  initially  zero. 

Some  special  problems  exist  with  certain  types  of  nonlinearities.  In  addition, 
the  control  system  equations  may  be  used  to  perform  certain  auxiliary  functions  such 
as  integration,  filtering,  or  engine  dynamics  modeling. 

If  gain  scheduling  is  used  on  any  transfer  function,  simply  set  the  element  gain  to 
unity  and  calculate  the  new  stage  input  k(t)  U(t)  in  the  previous  nonlinear  stage. 

Closed  loop  systems  containing  nonlinear  elements  present  a  special  problem. 


The  existence  of  nonlinear  elements  Inside  the  feedback  loop  prevents  reduction 
of  the  system  to  a  single  linear  element.  Since  both  u(t)  and  f(t)  must  be  known  to 
determine  y(t);  and  since  y(t)  must  be  known  to  determine  f(t),  y(t)  must  be  known  to 
determine  y(t)  without  such  a  reduction.  If  the  nonlinear  elements  are  schedules  or 
computational  nonlinearities,  they  change  the  system  characteristic  equation  and 
necessitate  the  calculation  of  the  state  space  matrices  in  real  time.  If,  however, 
they  are  switching  functions  there  are  only  three  possibilities; 
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1.  Switch  2  Is  open  (y(t)  = 

2.  Switch  2  is  closed,  switch  1  is  open  (v(t)  =  V^.^) 

3.  Both  switches  closed  (nonlinearities  ignorable). 

This  system  may  be  represented  by  three  separate  systems  over  the  next  cycle 
interval. 


NONL  IN  EAR  STAGE 

r"  VLIM - M 


The  succeeding  nonlinear  stage  may  compare  U(t)  ±  f(t)  and  to  determine 
if  Yj^(t)  or  Y2(t)  will  be  compared  to  thus  choosing  the  correct  value  for  Y(t) 

from  among  Y];^,  Yj^(t),  and  Y2(t). 

It  may  be  noted  that  very  simple  elements,  such  as  multiplication  by  a  constant, 
may  be  considered  either  a  zeroth  order  transfer  function  or  a  simple  computational 
nonlinearity. 

Because  of  the  time  required  to  execute  the  logic  of  the  linear  element  repre¬ 
sentation,  such  simple  elements  will  be  considered  most  efficiently  represented  as 
computational  nonlinearities. 

To  use  rotary  and  oscillatory  balance  data  for  departure  and  spin  modeling,  it 
is  necessary  to  separate  the  steady  and  unsteady  components  of  angular  rate.  In 
addition,  a  number  of  variables  describing  the  aircraft  flight  condition  must  be  inte- 
grared.  While  these  operations  are  not  explicitly  shown  as  part  of  the  control  system. 
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they  may  be  easily  represented  as  a  networic  of  linear  elements  just  as  are  the  linear 
stages,  i.  e. ,  the  subroutine  that  handles  the  linear  stages  may  be  called  from  the 
main  program  to  perform  these  functions. 

Reference  1  contains  several  F-14  control  systems  of  interest  which  have  been 
combined  here  to  form  control  systems  A  and  B.  ^stem  A  represents  the  current 
F-14  control  system  including  a  Grumman  designed  Aileron  Rudder  Interconnect  (ARI) 
intended  to  alleviate  wing  rock  and  roll  reversal  above  15-degrees  angle  of  attack. 

The  high  gain  of  this  system  has  been  found  to  create  lateral  pilot  induced  oscUlation 
(PIO)  tendencies  above  20-degree  angle  of  attack.  For  this  reason  the  F-14  is  cur¬ 
rently  being  operated  with  the  ARI  (Lateral  Stability  Augmentation  System  (SAS)  switch) 
off.  system  B  represents  a  design  effort  by  NASA  Langley  to  alleviate  the  PIO  ten¬ 
dencies  of  system  A  with  the  lateral  SAS  on,  and  has  recently  been  flight  tested  at 
NASA  Dry  den. 

The  following  diagrams  are  based  on  information  taken  from  references  1  and  2, 
but  they  have  been  rearranged  tc  facilitate  application  of  the  digital  control  system 
ret>resentation  presented  here.  The  engine  model  is  a  simple  first  order  filter  as 
used  in  reference  1.  No  attempt  is  made  here  to  model  the  Approach  Power  Com¬ 
pensation  system. 


25 


NADC-80220-60 


CONTROL  SYSTEM  VARIABLES 


SYMBOL 

SYSTEM  INPUTS; 

DEFINITION 

UNITS 

h 

Altitude 

ft. 

M 

Mach  No. 

— 

N 

z 

Normal  Load  Factor 

(a/g) 

N 

y 

Lateral  Load  Factor 

(V® 

-1 

p 

Roll  rate 

Sec 

Q 

Dynamic  pressure 

Lb/Ft^ 

-2 

q 

Pitch  acceleration 

Sec 

-1 

q 

Pitch  rate 

Sec 

-1 

r 

Yaw  rate 

Sec 

a 

Angle  of  attack 

Deg. 

Rudder  pedal  deflection 

In. 

*LAT 

Lateral  stick  deflection 

In. 

*LONG 

Longitudinal  stick  deflection 

In. 

^SASWT) 

Directional  SAS  switch  position 

— 

^SASWL 

Longitudinal  SAS  switch  position 

— 

^SASWLT 

Lateral  SAS  switch  position 

— 

*TRD 

Directional  trim  switch  position 

— 

®TRL 

Longitudinal  trim  switch  position 

— 

*TRLT 

Lateral  trim  switch  position 

— 

*THL 

Left  engine  throttle  position 

In. 

®THR 

Right  engine  throttle  position 

In. 

r 
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CONTROL  SYSTEM  VARIABLES  (Continued) 
SYMBOL  DEFINITION 

SYSTEM  INPUTS  (Cont'd. ) 


A  man 

A  man 
A 

^aut 

n 

SYSTEM  OUTPUTS 

^lat. 

^lor^ 

F  . 
ped 

^a  lim 

*e  lim 

lim 

^splim 

®tcl,  ®tcr. 


Position  of  manual  sweep  override  switch 

Manually  selected  sweep  angle 

Wing  sweep  angle 

Output  of  wing  sweep  schedule 

Total  angular  rate 

Total  force  applied  to  lateral  stick 
Total  force  applied  to  longitudinal  stick 
Total  force  applied  to  rudder  pedals 
Differential  tail  deflection  (limits  applied) 
Stabilator  deflection  (limits  applied) 

Rudder  deflection  (limits  applied) 

Spoiler  deflection  (limits  applied) 

Engine  thrust  levels 
All  other  variables  are  used  internally. 
Variables  E_  and  D„  represent  input  data 
associated  with  switching  function  and 
computational  nonlinearities  in  the  ith 
nonlinear  stage. 


UNITS 


Deg. 

Deg. 

Deg. 

Sec  ^ 

Lb. 

Lb. 

Lb. 

Deg. 

Deg. 

Deg. 

Deg. 

%  max  thrust 
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CONTROL  SYSTEM  DIAGRAMS 
SYSTEM 


FIGURE  la.  Control  System  A:  First  Nonlinear  Stage 
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TABLE  I. 

CONTROL  SYSTEM  A:  FIRST  LINEAR  STAGE 


FIRST  LINEAR  STAGE 


INPUTS; 


Ul 

U2 

U3 

U4 

U5 

Ue 

U7 

c 

00 

U9 

‘ILIM 

5  r2 

^LAT 

LIMSW 

_ 

^tre 

_ 

*tra 

_ 

^trr 

_ 

^MAN 

OUTPUTS: 


Vi 

V2 

V3 

V4 

vb 

ye 

y7 

Vs 

79 

5eq-| 

B 

^32 

B 

B 

^tre 

«tra 

*trr 

^MAN 

T ransfer  Matrix  entries: 


S  +  .5 


®22* 


-  8.56 

°33*sTr 


G44  = 


20. 

S  +  20. 


Gkk  - 


55  “  s  +  .5 


Gee* 


1 

s 


S 


G99  = 


1 

S 


Engine  dynamic  model 


“10 

“11 

“12 

“13 

^THL 

526 

®THr 

CD 

CM 

LU 

YlO 

yii 

V12 

Y13 

^THLI 

®THL2 

*THr1 

*THr2 

Gi010  =  G-,212  =  s^  *^11  11  =  ^313  13  “"rrs: 
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TABLE  II. 

CONTROL  SYSTEM  A. 
SECOND  LINEAR  STAGE 


INPUTS: 

“1 

“2 

“3 

“4 

“5 

“6 

“7 

3eq1  LIM 

«a3 

^32 

Sspc 

^34 

^’’B 

^36 

OUTPUTS; 

Vi 

- 1 

V2 

y3 

y4 

Vb 

ve 

y7 

to 

Sa^ 

5sp 

Ssp^ 

Srg 

Sr7 

Transfer  Function  Matrix: 

«  _  1.0114  {S+5)2  _  90 

(S+1.887)  (S+13.4)  ^22  "  ^66  ($+90) 

r  -  r  =90 
^33  -  ^77  ^ 

^  20  r  -  20 

^44  ~  S+20  '’SB  s 
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FIGURE  3.  Control  System  A;  Third  Nonlinear  Stage 
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TABLE  m, 

CONTROL  SYSTEM  A:  TfflRD  LINEAR  STAGE 


INPUTS: 

“1 

U2 

“3 

^eq2  LIM 

5rc 

^45 

OUTPUTS: 

Vi 

y2 

y3 

*eq3 

«r8 

■a 

Transfer  Matrix  Entries 
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TABLE  IV. 

CONTROL  SYSTEM  A:  FOURTH  LINEAR  STAGE 
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FIGURE  5.  Control  System  A:  Fifth  Nonlinear  Stage 
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.1 


CONTROL  SYSTEM  CONSTANTS 

TABLE  Va. 

CONTROL  SYSTEMS  A  AND  B 


FIRST  STAGE: 


Ell  =  -8727  sec"^  (50°/sec.) 

Ei7  =  2.356  sec'"'  (135“ /sec.) 

El  13  =  1  g- 

Ei2=  ISVsec. 

Ei8  =  19“ 

El  14  “  -8727  sec'^ 

Ei3  =  .5  in/sec. 

Eig  =  1  in. 

El  15=  1  in. 

Ei4  =  .19  in/sec. 

El  iQ  =  .55 

El  ig=  1  in. 

E-jg  =  .113  in/sec. 

El  11  =  1  in- 

El  17  =  .35  sec'^ 

Eie  =  4.363  sec'^  (250“ /sec.) 

El  12  =  1  in. 

El  18  =  44“ 

El  19  =  25“ 

El  20  =  "*5“ 

Dll  =  3lb./g. 

Dig  =  .15 

Di  11  =  -10“ /in. 

Di2  =  1  g- 

Di7  =  2. 

Di  12  =  9.15° /g. 

Di3  =  7.428  Ib-sec2 

Dig  =  -38.96° -sec. 

Di  13  =  22.5“ /g. 

0^4  =  2“ /in. 

Dig  =  -12.66“ /in. 

Di  14  = -57.3° -sec. 

Dig  =  .04 

Dl  10  =  '2°  /'"• 

SECOND  STAGE: 


E21  =  .2618  sec  "'  (15° /sec.) 

E23  =  10“ 

E25  =  21.5“ 

E22  = 

E24  =  50“ 

THIRD  STAGE: 


E31  =  .2618  sec'^  (15° /sec.) 

634=  12.5“ /sec.  /sec.) 

E32  =  .3667“ /sec.  (|§'’/sec.) 

^35  =  55“ 

E33  =  5° 

Egg  =  .7044“ /sec.  (^“/sec.) 

E37  =  19° 

NADC -80220-60 

TABLE  Vb. 

CONTROL  SYSTEM  CONSTANTS  (Contd. ) 


FOURTH  STAGE; 


E4t  =  .175  sec (lOVsec.) 
^42  =  3°  ^44  =  30“ 


E43  =  5.3“ /sec.  (^^“/sec.) 
D41  =  17.19“ -sec.  (.3“/ “/sec.) 


FIFTH  STAGE: 


E51  -  1.8“ /sec.  (||“/sec.) 

£53  =  35“ 

555  =  -33“ 

D5I 

E52  = 15“ 

554  =  10° 

E56  =  12° 

FOR  CONTROL  SYSTEM  B  CHANGE 


D^g  =  -6.33“ /in. 


Di8  =  -114.6“ -sec. 
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CONTROL  SYSTEM  SCHEDULES 

MACH  TRIM: 


S  LONG  - 


WING  SWEEP: 


A  ~  DEG. 


LONGITUDINAL  SPRING: 


^LONG  ~ 


LATERAL  SPRING: 


Flat  -  lbs 


^LONG  - 


Slat'' IN 


FIGURE  8a.  Control  Systems  A  &  B:  Schedules 
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FIGURE  8c.  Control  Systems  A  &  B;  Schedules  (Contd.) 
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GAIN  SCHEDULES 
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APPENDIX  A 

DIGITAL  MODELING  OF  CONTINUOUS 
LINEAR  ELEMENTS 


TABLE  OF  CONTENTS 
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DIGITAL  CONTROL  &  SIMULATION 

An  Important  difference  exists  between  digital  simulation  of  continuous  linear 
elements  and  the  design  of  digital  control  systems  even  though  they  both  employ  the 
mathematics  of  the  discrete  time  domain.  Digital  control  theory  is  concerned  with 
determining  the  sampled  or  continuous  response  of  a  continuous  system  subjected  to 
a  sampled  input. 

DIGITAL  CONTROL 


The  Z  transform  method  is  directly  applicable  and  may  be  used  to  obtain  exact 
solutions  because  the  actual  input  to  the  system  is  sampled.  This  result  is  achieved 
by  inverse  Laplace  transforming  the  transfer  function  to  obtain  the  impulse  response, 
and  then  Z  transforming  the  resultant  function  of  time  to  obtain  the  impulse  transfer 
function. 

G(z)  =  z[£  |G  (s)|] 

In  contrast,  digital  simulation  seeks  an  impulse  transfer  function  which  when 
subjected  to  the  sampled  input  yields  a  sampled  output  which  is  a  good  approximation 
to  the  sampled  version  of  the  continuous  output  signaL 

DIGITAL  SIMULATION 


The  equality  can  only  be  approximate  because  u(kTs)  represents  u(t)  only  at 
discrete  instants.  G'(z)  must  account  for  the  variations  in  u(t)  between  sampling 
instants.  Furthermore,  since  G(z)  represents  the  response  of  system  G(s)  to 
sampled  inputs  and  G*(z)  the  response  to  continuous  Inputs,  G(z)  and  G'(z)  are  not 
generally  equaL 
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For  computer  programming  purposes,  it  is  necessary  to  determine  G’(z)  and  to 
inverse  Z  transform  to  obtain  a  difference  equation  relating  the  output  at  any  sampling 
instant  to  the  previous  values  of  input  and  output.  To  fully  explain  how  this  is  done, 
a  brief  review  of  sampled  signals  and  Z  transforms  is  in  order. 

APPLICATION  OF  THE  Z  TRANSFORM 

A  sampled  signal  is  produced  by  an  ideal  switch  which  opens  and  closes  instan¬ 
taneously.  The  switch  usually  closes  at  regular  time  intervals,  t=kTg,  and  remains 
closed  for  a  veiy  short  period.  At.  This  period  is  sufficiently  small  that  the  sampled 
signal  may  be  considered  constant  during  the  sample  period. 

-uikTj) 


i__i - t 

t 

The  sampled  signal  is  therefore  a  series  of  square  pulses  of  finite  amplitude  and 

very  small  width.  The  time  interval  At  is  imprecisely  defined  but  constant.  As  such, 

it  represents  a  constant  scale  foctor  which  would  be  carried  along  through  any  linear 

equation  without  adding  much  to  the  understanding  of  the  system.  We  may  eliminate 

At  by  normalizing  the  square  pulse,  u  (t),  to  give 

s 

U*(t)  =  i  u^(t). 


u(t) 


Since  the  square  pulse  has  the  property: 
kTg  +  At 

kT/  Ug(t)  dt  =  u(kTg)  At. 

s 

The  normalized  pulse  has  the  property 

kTg  +  At 

y*  u*(t)dt  =  u(kTg). 


kT- 


From  this  equation  we  suspect  a  close  relationship  between  the  normalized  pulse  and 
the  unit  impulse  function,  whose  integral  is  unity,  and  which  is  comprised  of  two  step 
functions  whose  amplitudes  are  adjusted  as  At  —  0. 
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In  order  to  insure  u(t)  constant  over  the  Interval  At,  we  must  allow  At  to  become 
vanishingly  small  and  the  normalized  pulse 


Therefore,  the  normalized  output  of  the  sampler  is  a  series  of  impulses  whose  strengths 
are  equal  to  u(t)  at  each  sampling  instant. 

u*(t)  =  u(kT  )  5(t) 
s 

Laplace  transforming  gives 

“IcT  s 

jC  u*(t)  =  u(kT^)  e  s 

If  all  the  normalized  pulses  for  all  sampling  instants  are  Laplace  transfoimed  and 
summed,  and  if  the  substitution, 

T  s 
s 

z  =  e 

is  made,  the  expression  for  the  Z  transform  is  obtained 

9Q 

u(z)  =  Z  [u(t)]  =  u(kT^)Z 

k=o 

The  Z  transform  can  be  applied  to  any  function  of  time  and  may  be  viewed  as  a 
discretized  version  of  the  Laplace  transform 


The  Z  transform  is  a  special  case  of  the  Laplace  transform  where  the  function 
being  transformed  is  a  series  of  impulses;  therefore,  it  obeys  the  same  algebraeic 
rules.  It  also  has  some  unique  properties  of  its  own.  One  of  its  most  important 
properties  can  be  obtained  by  expanding  the  terms  of  its  definition  equation 
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u(T^) 

u(z)  =  u(o)  + -  +  .  .  . 


u  (nT  ) 

+  - —  + 

n 

z 


u((n+l)T  ) 
s 

ii+i 

z 


The  time  at  which  k=o,  t=o  is  an  arbitrary  choice.  Suppose  the  time  where  u(t)  =  u^Tg) 
were  chosen  as  the  time  for  k=o,  t  =  t+nTg.  The  Z  transform  of  this  function  expands 
to, 


u((n-H)T^) 

Z  [u(t+nT^)J  =  u(nT^)  + - - - + . 

whose  entries  differ  from  the  terms  of  u(z)  for  which  k>n  only  by  the  factor  Z^ 

z°  u(z)  -  Z  [u(t+nT  )]  =  z^  u(o)  +  z^  ^u(T  )  +  ...+  z  u((n-l)T  ). 

L  s  I  s  s 

THE  DIFFERENCE  EQUATION 


It  will  now  be  shown  how  this  property  allows  the  current  value  of  the  output  to  be 
expressed  in  terms  of  the  previous  values  of  input  and  output. 


Recalling  that  the  choice  for  t=o  is  arbitrary,  the  right  hand  side  of  the  preceding 
equation  may  be  made  to  vanish  if  the  system  was  at  some  previous  time  quiescent. 
This  may  be  accomplished  by  choosing  t«^ 


sufficiently  far  back  in  time  so  that  the 
first  n  samples  are  zero.  This  yields  the 
relation; 


z  u(z) 


=  Z 


u(t+nT 


u(t) 


2T,  3T,  4T,, 


More  importan^,  the  inverse  Z  transform  of  any  such  function  of  Z  is  simply  the 
value  u  at  the  n^sampling  instant  in  the  future. 


u(z)j  =  u(t+nT  ). 

s 


This  result  is  completely  general  and  applies  likewise  to  the  output  function,  y(t). 
In  general  y(z)  and  u(z)  are  both  infinite  series  but  may  be  represented  approximately 
by  a  finite  number  of  terms.  All  Impulse  transfer  functions  represent  the  Z  transform 
of  some  system's  impulse  response  which  does  not  change  with  time.  Therefore,  it  is 
always  possible  to  obtain  an  approximate  impulse  transfer  function  as  a  ratio  of  con¬ 
stant  coefficient  polynomials  in  Z. 


y(z> 

u(z) 


G(Z) 


Pjlz) 
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which  may  be  expanded  to 


P2oy  (2)+P2i(z)y(z)  +  P222^y(z)  +  •  •  •  +  ' 


PjO  u(z)  +  •  •  •  +  u(z) 


which  when  inverse  Z  transformed  gives 


m 


E 


P2i  y  (t+iT^)  = 


1=0 

or  alternately 


E 


1=0 


P  u  (t+iT  ) 
li  s 


n 


y  (t+mT 


E 


li 


u  (t+iT  ) 

S  X  J  P„  s 

2m  ,  ^ 

i=o  i=o 


m-1 

-E 


2i 


Simply  by  considering  the  value  of  y  (t+mT  )  to  be  the  current  value  of  the  output, 
we  can  make  the  substitution 


t  =  t  -  mT 


to  obtain 


y(t)  = 


E 


1=0 


li 


u  (t-(m-i)T  ) 
^2m  ® 


m-1 


E 


1=0 


2i 


y  (t-(m-l)T  ) 
^2m  ® 


We  have  the  current  value  of  the  output  in  terms  of  the  previous  values  of  input  and 
output.  However,  a  method  for  determining  the  approximate  impulse  transfer  function 
G’(Z)  in  terms  of  the  ratio  of  constant  coefficient  polynomials  in  Z  must  still  be 
developed. 

CONVERSION  OF  TRANSFER  FUNCTIONS  TO  DIFFERENCE  EQUATIONS 

A  linear  element  is  customarily  defined  by  a  transfer  function  that  is  usually 
expressed  as  the  ratio  of  factored  polynomials  in  S. 
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m  >  h 


If  the  transfer  function  pol5momials  are  expanded  and  divided  by  S^,  this  equation 
may  be  stated  in  terms  of  the  forms 


a  y(s)  b  u(s) 

—  and  — 


Inverse  Laplace  transforming  these  terms  is  equivalent  to  integrating  a^  y(s),  j  times 
or  bj[  u(s),  i  times.  If  the  integral  can  be  Z  transformed  as  a  pol3momi^  in  Z  times 
u(z)  or  y(z),  a  representation  of  G'(Z)  of  the  required  form  will  be  obtained. 

The  sampled  signal  y(z)  represents  y(t)  only  at  discrete  instants  t  =  kTg.  In 
order  to  integrate  y(t)  some  approximation  to  y(t)  between  sampling  instants  must  be 
assumed. 


rectangular 


trapezoidal 


To  obtain  the  value  of  the  integral  in  poljmomial  form,  y(t)  must  be  expressed 
(approximately)  as  a  polynomial. 
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Where  the  constants,  g^,  are  determined  by  satisfying  the  constraint  equations 
6 

h  k  >  j  >  k  -  C 

i=o 

(k  =  latest  sample  no. ) 

The  integral  might  be  determined  by  considering  the  polynomial  to  represent  y(t) 
over  either  of  two  intervals; 

easel;  kT  ^t>(k-l)T 
s  s 

case  2:  (k+l)T  >  t  >  kT 
s  s 


Case  1  would  be  expected  to  be  more  accurate  for  ?>0,  because  y(t)  at  both  ends  of  the 
interval  are  known.  For  <!=0  there  seems  to  be  no  particular  advantage  to  either  case. 

The  choice  does  make  a  difference  in  the  final  form  of  the  Z  transform.  For  now 

we  will  refer  to  the  integration  limits  as  kT  and  (k-l)T  ,  so  that 

s  s 


U  +1  .  ,  T+1 

^  -  (k-1) 


over  all  j  intervals  from  1  to  k 


i=o 
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Z  transforming,  and  realizing  that  /  y(t)  dt  =  0 


/ 

•/  n 


Case  1 


/ 


kT 


y(t)  dt 


k  e 


SEE  4 


ii  T 

i+l  s 


.i+1  i+l  -k 

j  -  0-1)  z 


k=l  j=l  i=o 
For  case  2,  the  limits  are  changed  so  that. 


/. 

kT 


(k+l)T^ 


y(t)  dt  = 


■  Z-f  i+l 


T 


i+l 


(k+1)^^^  -  k^^^ 


S  1=0 

But  there  are  now  k  +1  intervals 


I 


(k+l)T  k  ft 

j=o  i=o 


^ij  „  i+l 
i+l  s 


,  i+l  ,i+l 

0+1)  -  J 


Also  the  value  of  the  integral  is  not  zero  at  k=o,  since  /  y(t)dt  ^  0 


/' 


Case  2 


r  (k+l)T^ 

'  E  E  E 

z 

J  y(t)  dt 

**  0 

o 

II 

o 

II 

1? 

-k 


Although  not  readily  apparent,  this  infinite  series  is  equivalent  to  a  geometric 
progression  expressible  as  a  finite  polynomial  in  Z.  Since  specific  expressions  are 
needed,  it  is  best  to  demonstrate  this  fact  with  an  example. 
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e=0.  y(t)  =  g  .  =ya'T  ) 

v/J  O 


CASE  1 


Consider  the  simplest  rectangular  or  zeroth  order  approximation  to  y(t). 

rectangular  approximations 

The  Z  transform  of  the  case  1  approximation  y 

gives 

kT 


L^fTTTr 


/, 


y(t)  dt 


k=l  j=l 


y(jTg)Ts  z 


-k 


y(t) 


Whereas  case  2  gives 
.(k+l)T 


/ 

**  r\ 


y(t)dt 


y(JT^)  z 


-k 


k=o  i=o 


Expanding  the  summations 


/•kT 

/  y(t)  dt 

Jo 

y(T)  ^  y(T)  +  y(2T)  ^  _ 

z  2 

z 


/: 


(k+l)T 


y(t)  dt 


=  T, 


y(o)^y<°>ty0  +  ... 


CASE  2 


We  see  the  difference  between  the  Z  transform  in  cases  1  and  2  is 

Tgy(o)  ji  ^  +  .  .  .| 

The  bracketed  term  is  a  geometric  series  which  we  can  see  has  a  finite  expression 
if  we  allow 
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For  an  initially  undisturbed  system  y(o)  =  o,  and  the  Z  transform  of  the  integral 
is  independent  of  the  integration  limits. 


[/y(t)dt]  = 


TsZ 

(z-1) 

Similarly  for  case  1  with  8=1 


y(z)  «=0. 


/. 


kT 


y(t)dt 


■’.ri; 

K=1  j=l 


Also;  y  OT  )  =  ) 

s  oj  i]  s 


ya-i)V  =  goj^Sij<(^-^>V 


so  that. 


g„j  •  J  y(a-l)T^I  -  (l-l)  y  OT^) 


Therefore 


I 


kT 


y(t)  dt 


=  T, 


r . 


k=l  j=l 


-  14^43^)  +  yOT  ) 

S  5  S 


-jy4=A^  +  jyliX^-  1/2  y(jT  )  >  Z 

So  S 

+  1/2  y(O-l)Ts) 

V.  J 


-k 


irZE  ( 

k=l  j=l 


y  ((j-i)Tg)  +y(jT^)  J 


l-k 


Expanding 


r 


y(o)  +y(T  )  +  y(o)  +  2y(T  )  +  y(2T  ) 
s  _ s  s 


+  y(o)  +  2y(T  )  +  2y(2T  )  +  y(3TJ 

S  o  o 


+  .  •  •  • 
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Up  m-1  p 

y(t)  =y;  u(t-(ni-i)T  )  -  ^  y(t-(in-l)T  ). 

"  ^2m  ®  "  ^2m  ® 

1=0  1=0 

USE  OF  POLYNOMIAL  EXPANSION  AND  THE  TUSTIN  TRANSFORMATION 

To  relate  the  coefficients  of  this  equation  to  the  transfer  function  roots  an(^  to  the 
order  of  the  approximation  to  y(t)  and  u(t),  we  need  only  have  a  general  method  of 
expanding  polynomials.  Such  a  method  can  be  found  by  expressing  a  polynomial  after 
i  factors  have  been  multiplied  out  as; 

2  ^lii  ®00  " 

i=o 
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All  additional  factors  may  be  broken  into  first  or  second  order  factors .  If  the 
next  term  is  first  order, 


Pj^l(s)  =  (s  +  1/r)  Pj,  (s)  = 


8+1 

L 

i=0 


o  i  ^ 
(«+l)  ® 


8+1 


8+1 


a,l  =  0 
)  if  i  >  8 
or  i  <  0 


If  the  next  term  is  second  order 


Pj+2(s)  =  (3^+2  +  P^(s) 


8+2 


P,+2(S)  = 


2  1 
7  (a.  u  +  2  (w  a  +  a 

X  j  '  gi  n  "  u  D/,‘_o\/  ° 


n  8  (i-1)  8(i-2)' 


i=o 


If  we  make  the  substitution 


s 


into  the  n  order  polynomial  is  S,  we  obtain  an  8n  order  polynomial  in  Z,  and  we 
must  keep  track  of  the  last  8n  values  of  input  and  output  to  calculate  y{t).  There  is 
therefore  a  tradeoff  between  complexily  and  accuracy  in  selecting  a  value  for  8 .  The 
substitution  for  8=1  has  been  used  extensively  as  a  good  compromise  value  and  is 
usually  called  the  Tustin  transformation.  * 


s 


2  (z-1) 

Tg 


If  a  transfer  functions  n  order  numerator  and  m  order  denominator  are  ex¬ 
panded  into  polynomials  of  s,  the  Tustin  transform  may  be  applied  to  give 


*See  reference  4 
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n 

u(z)k  2  b. 
1=0 


i 


where  k  is  the  transfer  function  gain. 


Multiplying  both  sides  of  this  equation  by 


(z+1) 


m 


gives: 


„i  m-i  ,  i  , 
a.  2  T  (z-1)  (z+1) 
1  s 


n 

k'U(z)^^J  b^  2^  ^  (z-1)^  (z+l)°^  ^ 

i=o 


th 

The  result  is  that  both  sides  of  the  equation  contain  m  order  polynomials  in  z.  If  the 
expansion  of  the  i*  set  of  z  factors  is  defined  as 

m 

V”'  1  i  .  m-i 

y  h..  z-*  =  (z-l)  (Z+1) 

j=o 

h^j  for  all  values  of  m  >  i  >  o  may  be  calculated  using  the  factor  expansion  process 
developed  on  page  A- 14. 


This  then  makes  the  expanded  version  of  the  impulse  transfer  function  (G'(z)) 


m  m 


i=o  i=o 


n  m 


i=o  j=o 


A-15 


NADC-80220-60 


Making  the  coefficients  of  the  difference  equation 


y(t) 


-JL  y(t-(m-j)T  ) 
^2m  ® 


eqxial  to 


n 


i=o 


Thus  the  difference  equation  representing  the  inverse  Z  transform  of  G'  (z)  may 
be  obtained  directly  from  the  transfer  function  G(s). 

In  order  to  calculate  y(t)  directly  using  this  formulation,  we  require  the  storage 
of  2m-l  coefficients  and  a  like  number  of  past  input  and  output  values  for  each  element. 
Significant  efficiencies  in  representing  single  and  multiple  element  systems  can  be 
achieved  using  the  state  space  notation  of  appendix  B.  This  method  minimizes  the 
variables  and  algebraeic  manipulations  required  to  calculate  all  system  outputs  while 
also  presenting  a  very  flexible  means  of  representing  the  system  element  structure. 

In  order  to  demonstrate  this  formulation  two  simple  examples  are  now  presented. 
Example  1;  First  order  filter 

Using  the  expansion  procedure 
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IS 
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“lo  ■  *  \a  " 


=  0 


=‘-  -“-*>15?””^ 


11  00 
Also  expanding 


(Z-1)'  (Z+l)^"' 

m 

=Z  Nr 

j=o 

h  =1 

K,  =  1 

00 

01 

^0  =  -^ 

II 

Evaluating  the  difference  equation  coefficients 


P  =  b  T  h  =  T 
10  00  S  00  S 


P  =  b  T  h  =  T 
11  00  S  01  S 


'’20**10’^S''00*^1'''10*<'’V'> 


P  =a  Th  +a.2h  =  (aT  +2) 

21  10  S  01  11  11  '  S  ' 


making  the  difference  equation. 

k 


k  ,  1 


Example  2;  First  over  second  order  system 

_ ,  k  ( s  +  o)  _ 

G(s)  =  —2 — ^ - j-  n=l,  m=2 

(s  -*■  2  f  3  +  ) 


Expanding  the  numerator  polynomial  as  with  example  I's  denominator 
^o  ■  ^l  =  *• 
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Expanding  the  denominator 


a  =  1 
00 


^  '  |a  ,  cj  ^  +  2  fcj  a  +  a  | 
2  (  oi  n  n  o  i-l  oi-2i 


“20  ■  “OO  ^  *  >6^2  ■ 


=  ^  +  2  f  cj  a  +  =  2  icj 

n  n  00  -^-1  n 


2 

+2  +  a. 1 

22  x^2  n  'n  01  00 


Expanding 


(z-l)‘  (z+l)“"^  =  h.^  Z^ 


2>i>0 


for  i  =  o 


(Z+1)  (Z+1)  =  >  h^.  Z^ 


using  the  notation 


eh  =  h  after  e  factors  expanded 
i]  i] 


“oo  =  ' 

“■oo '  "  no  *  * 

^”01  “  "'■oo  "  <3^  *  ‘ 


1 

J 
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=  2h 

00 


00 


■  %  *  <^“00  *  " 


\l  ^0  "  ''“oi  *  " 

"oa  ■  '“02  '  %  * * i>  -  ‘ 

which  can  be  verified  by 
2  2 

(Z+1)  =2  +  2z  +  1 

Without  bothering  to  demonstrate  the  expansion  formula  for  i=l  and  i=2  we  have 


[  =  -1 

h  =0 

h 

=  1 

10 

11 

12 

=  1 

h  =  -2 

h„„ 

=  1 

20 

21 

22 

The  difference  equation  coefficients  are; 


^10  ~  ^10  ^00  ^11  ^  ^10 


=  0X3  -  2X3  =  T3  (0X3-2) 


P,,  =  X^  h„  +  b  2  X  h  , 
11  10  S  01  11  S  11 


=  2aX„ 


P  =b  X  h  +b  2Xh 
12  10  S  02  11  S  12 


=  0X3  +  2  Xg  =  Xg  {oTg+2) 


P  =a  X  h  +a  2Xh  +a  4h 
20  20  S  00  21  S  10  22  20 


2  2 

=  CJ  -  4  f  4 

n  S  S  n 
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^21  ^20  ^21  ^  ^11  ^  ^22  ^21 


=  2  (oj  ^  T^  -  4) 
n 


^22  ■■  ^20  '^S  ^02  ^  ^21  ^  ^S  ^12  ^22  ^  ^22 


4?  T 

2^2  s 

=  u)  T_  +  0)  +4. 

n  S  n 


Making  the  difference  equation 

y(t)  = 


k  Tg  (oTg-2) 


4  f  T 

2  2  *  S 

to  T„  +  — C3 -  +  4 

n  S  n 


u  (t-2T  ) 
'  S 


2k  aT, 


4  f  T 

2  _  2  ^  ‘ 

w  T  +  -  +  4 

n  S 


u  (t-Tg) 


k  (ciTg+2) 


^  4  f  T 

2  2  S 

u>  +  4 

n  S  n 


u(t) 


^  tJ  +  4  ^  TVw  +  4 
n  S  S  n 


y  (t-2Tg) 


y  (t-Tg) 
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APPENDIX  B 

APPLICATION  OF  STATE  SPACE 
NOTATION 
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THE  STATE  SPACE  CONCEPT 


Any  time  invariant  linear  system  may  be  described  in  the  continuous  time  domain 
by  a  linear  ordinary  differential  equation  of  the  form. 


^4. 

dt 


,m-l  ,  ,n 

d  V  dv  .  d  u 

a, - f+...+a  ,-^+a  y  =  b - +*->+bu  n^m 

m-1  dt  m  o  ,.n  n 

dt 


^  dt“-' 


or  in  the  discrete  time  domain  by  a  difference  equation  of  the  form. 

P  P  kP  kP 

yW  +  f  y  +•  •  •  +  y  (t-mTs)  =  U(t)  +•  •  •  u(t-mTs) 


2m 


■  2m 


■  2m 


'  2m 


The  similarity  of  these  forms  suggests  that  either  equation  may  be  considered  a 
linear  combination  of  2m  dependent  variables.  Those  in  the  continuous  domain  are 
related  by  differentiation,  those  in  the  discrete  domain  by  a  time  shift.  The  analogy 
is  complete  in  the  Laplace  and  Z  domains  since 


^m-1  dt 


=  a  s  y(s)  and  z 
m-l 


^  y  (t-(m-i)Tg) 
2m 


_ 2i  i-m  ,  ^ 

^ —  z  y  (z) 

^2m 


Thus  realizing  that  we  may  interchange 


i'  i-m 
s  ■'^^z  .  a. 


■  2  m-i 


kP 


1  m-i 


,  etc. 


2m  2m 

we  may  adopt  a  notation  convention  for  these  two  equations  such  that  for  any  function  g; 


=  if  or  gj  =  g  (t  -  (m-j)  T  ). 
dt^ 


j  =  m-i 


The  differential  and  difference  equations  now  have  a  common  form  which  may  be 
expressed  in  either  of  two  useful  forms.  The  coefficients  of  the  differential  equations 
are  used  for  simplicity. 


.  m  ,  m 

(a  y-b  u)  +  (a  y»-b  u')+...  +  (y  -b  u  )  = 
m  m  m-l  m-l  o 


Eq.  .1 


B-1 
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or 


Output  State  Input 


£q.  .2 


The  output,  y,  is  dependent  only  on  the  input,  u,  and  the  first  m  higher  orders  of 
y  and  u.  If  we  are  only  interested  in  the  output  value,  we  may  set  all  orders  higher 
than  to  any  value,  especially  zero,  without  affecting  the  output  calculation.  In 
doing  so,  we  cause  equations  1  and  2  to  imply  a  whole  system  of  m+1  equations 
obtained  by  differentiating  or  shifting  equation  1  m  times. 

Equation  Set  A 


m 


b  u)  +  (a  y^-  b 
m  m-1  m-I 


1  m  ,  m  ^ 

u)+«..  +  (y  -bu  )=0 


my 


-  b  u^)  +  (a  ,  y^  -  b  u^)  +  •  •  •  +  (a,y°^  -  b  u™)  = 
m  m-1  m-1  1  1 


m-1  m-1^  ,  ,  m  ,  „ 

m  m  m-1  m-1 


m  ,  m 
ay  -  b  u  =0 
m  m 

In  the  case  where  set  A  represents  difference  equations,  it  is  clear  that  the 
various  orders  of  y  and  u  represent  the  past  history  of  the  input  and  output.  The  same 
is  true  for  the  differential  equation  interpretation  of  set  A  since  the  derivatives  of  y 
and  u  represent  those  time  histories  by  virtue  of  their  Taylor's  series  expansions. 

Because  an  initially  quiescent  system  requires  some  type  of  input  to  disturb  the 
system  and  produce  an  output,  output  may  intuitively  be  considered  a  direct  conse¬ 
quence  of  input.  But  the  system  equation  indicates  the  output  also  depends  on  the  past 
history  of  both  input  and  output,  by  virtue  of  the  higher  order  terms  in  y  and  u. 

This  intuitive  notion  leads  to  the  concept  of  system  state  as  a  combination  of  all 
factors,  distinct  from  the  input  and  associated  with  past  events,  which  contribute  to 
the  output.  The  state  may  further  be  defined  by  a  set  of  state  variables  whose  values, 
along  with  that  of  the  input,  determine  the  output.  It  is  clear  from  equation  2  that  state 
variables,  x.,  may  be  defined  so  that 

y  =  I  C.x  +Du 
1  i 
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which  is  most  often  written  in  matrix  notation  as 
y  =  C  X  +  Du 

where  C  is  a  row  vector,  D  is  a  constant,  and  x  a  column  vector  (the  state  vector) 
whose  components  are  the  state  variables.  What  is  neither  clear,  nor  in  fact  unique, 
is  how  to  associate  the  state  variables  with  the  terms  of  equation  1. 

A  logical  first  step  in  establishing  such  an  association  is  to  determine  how  many 
state  variables  are  required.  Equation  set  A  represents  all  the  equations  necessary 
to  describe  the  system.  Since  there  are  m+1  equations,  m+1  variables  are  required 
to  make  the  set  determinate.  That  is,  set  A  may  imply  the  value  of  an  m+2'^‘^  variable 
if  it  can  be  restated  in  terms  of  m+1  other  variables  whose  values  are  toown.  Since 
y  is  the  variable  whose  value  is  sought,  it  may  be  considered  the  m+2'^°  variable 
whose  value,  according  to  the  state  space  form  of  equation  2,  is  determinable  directly 
from  u  and  the  state  variables.  Therefore,  there  must  be  a  minimum  of  m  state 
variables,  and  any  additional  state  variables  are  unnecessary  and  may  be  eliminated. 

Consequently,  if  they  are  properly  defined,  there  wiU  be  m  state  variables  which 
are  linearly  independent  since  none  may  be  eliminated.  Therefore,  the  state  space 
notation  is  more  efficient  than  a  straightforward  application  of  the  difference  equation 
because  it  reduces  the  number  of  variables  (and  their  associated  manipulations)  from 
2(m+l)  to  m+2. 

The  state  space  form  of  equation  1, 
y  =  Cx  +Du 

may  be  considered  a  partial  definition  of  the  state  variables  which  specifies  the 
relationship 


But  one  such  equation  is  not  enough  to  define  the  constant  D  and  the  m  state  variables 
which  do  not  appear  in  equation  set  A.  An  additional  m  equations  are  required  which 
must  not  only  insure  the  independence  of  the  state  variables  but  must  be  consistent 
with  set  A  in  order  to  represent  the  systems.  That  is,  x  and  D  may  be  defined  by  any 
set  of  m+1  equations  that  are  linearly  independent  and  related  to  each  of  the  variables 
of  set  A  by  some  linear  combination. 

This  principle  allows  many  possible  formulations  of  the  additional  definition 
equations.  Highly  arbitrary  choices  may  be  involved  in  the  selection  of  any  particular 
form.  Simplicity  then  becomes  the  guiding  principle  in  such  choices. 
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If  we  retain  the  form  of 
y  =  Cx  +Du 

for  the  additional  m  equations  we  have 


fj^  =  +  BjU 


f  =  A  X  +•  B  u 
mm  m 

or  in  matrix  form, 

F  =  A  X  +  Bu 


The  independence  of  the  state  variables  may  be  assured  by  requiring 
det  A  /  0 

This  follows  from  the  property  of  determinants  which  allows  any  constant  multiple  of 
one  row  to  be  added  to  another  without  changing  the  value  of  the  determinant.  If  any 
state  variable  is  dependent,  it  may  be  expressed  as  some  linear  combination  of  other 
state  variables.  It  then  would  be  possible  to  manipulate  the  rows  of  A  in  this  manner 
until  either  one  row  or  one  column  is  all  zeroes.  Again  by  the  properties  of  deter¬ 
minants,  det  A  =  0  under  these  conditions. 

To  insure  the  definition  equations  are  consistent  with  set  A,  we  require  each 
variable  of  set  A  be  expressible  as  a  linear  combination  of  the  variables  requiring 
definition  (i.  e. ,  Xj  and  D).  If  v  represents  the  variables  of  set  A,  w  the  variables  of 
the  definition  set,  and  Q  the  matrix  relating  the  two. 


V  =  Qw 

Furthermore,  if  Pv  represents  equation  set  A,  each  equation  of  set  A  may  be  directly 
equated  with  a  linear  combination  of  Xj  and  D. 
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A  very  straightforward  way  of  iosurit^  consistency  is  to  impose  a  direct  corre¬ 
spondence  between  the  equations  of  set  A,  and  the  lines  of 


y  =  Cx  +  Du 
F  =  Ax  +  Bu 


definition  set. 


In  equation  set  A,  each  line  is  obtained  by  incrementing  the  equation  order  and 
dropping  out  terms  higher  than  m^^  order.  The  same  must  be  true  of  the  definition  set. 


Equation  Set  B 


1  ^-1  1  .  .  ~ 

-y  +Cx  +  Du  =-f^+A^x  +  B^u 

-  f^,  +  A,  x^  +  B,u^  =  -  f  +  A„  X  +  B  u 
11  1  2  2  2 


f^  ,  +  A  ,x^  +  B  ,  u^ 
m-1  m-1  m-1 


f  +  A  X  +  B  u 
m  m  m 


Noting  the  similarity  between  the  terms  on  the  night  above,  and  the  original  definition 
set,  we  can  make  a  judicious  choice  for  f^. 


=  y^  -  D  u^ 

h  "  ^1^  -  u^  =  y^  -  Du^  -  B^u^ 


f  =f  ,u=y  -Du  -  B,u  -  ...  -  B  u 

m  m-1  m-1  1  m-1 


which  simplifies  the  form  of  equation  set  B  to 


C 

A, 


m-1 


=  A  X  +  Bu 
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Further  choosing 


"c 

■  1 

A1 

=  I  = 

1  0 

• 

• 

0 

A  , 

1 

m-1 

Simplifies  the  form  for  the  definition  set 
y  =  C  X  +  Du 

=  Ax  +  Bu. 

If  this  system  represents  a  continuous  system  (differential  equation),  we  have  the 
familiar  forms: 

y  (t)  =  C  X  (t)  +  D  u(t) 

x(t)  =  Ax(t)  +  Bu(t) 

If  the  system  is  discretized  (difference  equation)  we  have 
y  (t  -  mT  )  =  C  X  (t-mT  )  +  D  u  (t-mT  ) 

S  S  5 

x(t-  (m-1)  T  )  =  Ax  (t-mT  )  +B  u  (t-mT  ) 
s  s  s 

Shifting  t  =  0  forward  by  m  sampling  intervals  yields  the  more  customary  forms. 


y  (t)  =  C  X  (t)  +  D  u  (t) 


X  (t  +  T  )  =  A  X  (t)  +  B  u  (t) 
s 


Knowing  these  forms  the  definition  equations  may  be  rewritten. 


x^  =  y  -  D  u 

x^  =  B^u  =  y^  -  D  u^  -  Bj^u 


1  m-1  ^  m-1  „  m-2 

X  =  X  -  B  u  =  y  -  Du  -  B  u 

m  ra-1  m-1  l 


-  B  u 
m-1 


B-6 


We  now  have  arrived  at  the  most  common  state  space  formulation.  * 
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differently,  the  form  of  B  could  also  be  simplified  to  express  directly  the  coefficients 
of  the  orders  of  u  above.  That  is. 


®i  = 


We  may  maintain  the  correspondence  of  the  rows  of  the  definition  set  with  the 
equations  of  set  A,  while  obtaining  the  desired  form  for  B  by  collecting  all  higher 
order  terms  together  as  a  single  first  order  state  variable. 


a  X 
m  1 


b  )  u 
m 


+  X 


=  0 


where 


1  1  U  V  1  ’"'1 

{,=a  ,x+(a  ,)u+...+x, 

k  m-1  1  m-1  0  m-1  1 


decrementing  the  superscripts 


\  =  v-1  ^  <"m-i  %  -  - 

^  k-1 

We  now  may  continue  defining  the  state  variables  in  this  manner  until  finally  we  get 
\  -  (m-2)  *  “l-’'!  ^  *  "'l 

We  now  have  a  system  of  m  + 1  equations  each  of  which  was  obtained  by  decrement¬ 
ing  the  superscripts  of  the  previous  one.  Therefore,  to  maintain  the  ordered  corre¬ 
spondence  with  the  lines  of  set  A  which  were  obtained  by  incrementing  the  superscripts, 
we  must  reverse  the  order  of  these  equations  setting  k  =  m  so  that 


X  =-a  X  +x  +fb  -ab)u 

1  11  2  ^1  10' 


^  2  "  ■  ®2  '‘l  ^  ‘'^3  ^  ^^^2  "  ^2%^  “ 


X  ,=-a  ,x,  +  x  +(b  ,-a  ,  b„)  u 

m-1  m-1  1  m  m-1  m-1  0 


=-a  x+(b  -a  b  )  u 
“  ml  m  m  0 


m 
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Since  the  latest  form  of  equation  1  on  page  B-8,  is  predicated  on 


y  = 


+  b  u 
0 


The  values  of  C  =  (1,  0. . .  0),  D  =  bQ  are  unchanged  while  the  A  matrix  is  altered 
slightly  and  the  B  vector  is  simplifi^. 


— 

— 

— 

— 

1 

0 

0 

A  = 

-^2 

0 

1 

...  0 

B  = 

^2 

-  a  b 

2  0 

. 

•1 

• 

-a 

m 

0 

0 

...  0 

b 

m 

m  0 

- 

- 

The  substitutions  of  page  B-1  may  be  made  for  discrete  systems. 


MULTIPLE  INPUT  AND  OUTPUT  SYSTEMS 


Until  now,  applications  of  state  space  notation  only  to  single  ioput,  single  output 
systems  have  been  considered.  The  concept  can  now  be  extended  to  a  network  of  linear 
elements.  First  consider  a  system  of  i  uncoupled  linear  elements  each  having  its  own 
set  of  state  space  equations  of  the  form, 


y. 

1 


X 

li 


+  b 


Oi 


u. 

1 


x/  =  A  X.  +  B.  u. 
1  i  1  11 


Combining  all  8  such  sets  into  a  single  pair  of  matrix  equations  gives 


y  =  X, 


X 


1 


* 


=  A 


X  +  B  u 
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where  the  dimensions  of  the  equations  have  been  expanded  so  that 


^11 

^12  '  ’  ’ 

^2 

A 

X  = 

^21 

X22  .  .  . 

• 

m 

1 

1 

1 

1 

^ml 

^m2  •  •  • 

^m2 

v 

< -  2  - ^ 


^2  •  •  • 

^2 

B  * 

®22 

®22 

1 

_®ml 

®m2  •  • 

•®m2 

^ _ 

.  P  _ 

m  =  order  of  highest  order  element. 

The  columns  of  X  are  the  state  vectors  of  the  2  elements  and  the  columns  of  B 
are  the  corresponding  input  vectors.  The  matrix  A*  always  has  a  constant  form, 
depending  on  the  expressions  chosen  for  the  B  matrix  entries,  of  either 


'^li 

1 

0  .  .  . 

O'" 

0 

1 

0 

.  .  0  ■ 

A*  = 

'hi 

0 

1  .  .  , 

0 

or 

0 

0 

1 

0 

• 

• 

1 

■  .  1 

-a  . 
mi 

0 

0  .  . 

0 

_ 

• 

• 

:  : 
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“multiplication  A  X  is  carried  out  in  the  usual  row  by  column  fashion  except 
that  ajj  is  selected  from  the  i®“  row  of  the  matrix 


^11 

CM 

^ml 

A 

^12 

^22  ■  •  • 

^m2 

• 

• 

• 

^1« 

"2  •  •  • 

a 

mv 

< - —  m 


when  multiplying  any  row  of  A*  by  the  i^^  column  of  X. 

Consider  now  the  case  in  which  the  «  linear  elements  are  coupled  by  pure  gains. 
The  inputs  to  the  elements,  Q*,  may  be  expressed  as  8  linear  combinations  of  k  sys¬ 
tem  inputs,  a. 


SIMPLY  COUPLED  INPUT  AND  OUTPUT 


That  is 
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The  output  may  be  similarly  coupled  giving  (for  j  system  outputs,  y) 


^11 

*^21  • 

r 

•  • 

.  .  C 

12 

22 

j 

• 

• 

• 

• 

• 

• 

. 

. 

. 

.  -  C  . 

y 

_  1] 

^ - 

2] 

D 

Substituting  into  the  multiple  element  state  space  equations 
7  =  Cx^  +  D  u 

=  A*X  +  B  u 

where 


D  =  C  b^  E 
B  =  B*  E 


and  b’  is  the  B  matirx  on  pg.  B-11. 

For  more  complex  types  of  coupling,  block  diagram  algebra  may  be  used  to 
reduce  the  network  to  this  simple  form.  Some  simple  examples  are; 
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OR: 


=> 


Gil 


It  is  clear  that  each  coupling  element,  including  pure  gains,  may  be  thought  of  as 
creating  a  separate  entry  in  the  transfer  matrix. 


'll 

S2-  • 

•  ^Ik 

G  .  . 

.  G  , 

'21 

22 

2k 

r 

G  .  .  . 

.  G . 

ji 

12 

jk 

This  means  that  any  network  may  be  reduced  so  that  entries  of  matrices  C  and  E 
are  either  0  or  1.  In  the  case  of  matrix  E,  only  one  entry  in  each  row  may  be  nonzero. 
Furthermore,  which  entries  are  nonzero  may  be  keyed  to  the  subscripts  of  the  trans¬ 
fer  function;  i.  e.  its  location  in  the  transfer  matrix. 


If  Go  =  G,,  , 
jk 


=  1 
k6 

otherwise  E. 

le 

=  0 

II 

Si 

=  0 

All  that  is  required  to  specify  a  highly  general  system  by  the  method  presented 
here  is  to  specify  the  number  of  elements,  their  location  in  the  transfer  matrix,  and 
the  gains,  poles,  and  zeroes  of  each  element.  The  state  spare  matrices  A*,  P,  C  and 
D  may  be  calculated  directly. 
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